function [g_num,perc_unk,trans_p] = bsoid_mt(data,fps,pix_cm,smth_hstry,smth_futr,min_trans) % data{m}, pixel/cm
%BSOID_MT     Classifies 3D pose estimation (DeepLabCut) output for rodent open field behavior using manual thresholding criteria.
%   
%   [G_LABEL,G_NUM,PERC_UNK] = BSOID_MT(DATA,PIX_CM) outputs classified behaviors based on DeepLabCut analysis
%
%   INPUTS:
%   DATA    A cell-array containing either one or multiple rodents' estimated positions of the 6-body parts (snout, 4 paws, and base of tail)
%           over time recorded from the bottom-up perspective. Rows represent chronological frame numbers. Columns represent x and y coordinates of 
%           individual body parts as follows:
%           Columns 1 & 2 tracks snout; columns 3 to 6 tracks the two front paws; columns 7 to 10 tracks the two hind paws;
%           columns 11 & 12 tracks the base of the tail. Tested on tracking data generated by DeepLabCut 2.0.
%   FPS    Rounded frame rate, can use VideoReader/ffmpeg(linux command) to automatically detect the input video fps.
%   PIX_CM    Pixel-to-cm conversion to adapt our manual thresholds to the user's set-up.
%   SMTH_HSTRY    BOXCAR smoothing using number of frames from before. Default ~40ms before.
%   SMTH_FUTR    BOXCAR smoothing using number of frames from after. Default ~40ms after.
%   MIN_TRANS    Minimum number of frames for transition to occur. Default ~100ms. 
%
%   G_NUM    Classified behavioral state, numbered.
%   PERC_UNK    Percentage of total number of frames that are undefined.
%   TRANS_P     Transition probability matrix with Markov assumption (rows as current, columns as next).
%
%   Created by Alexander Hsu, Date: 051519
%   Contact ahsu2@andrew.cmu.edu

    if nargin < 3
        error('Please input dataset, frame rate, & pixel-to-cm conversion!')
    end
    if nargin < 4
        smth_hstry = round(0.05/(1/fps));
        smth_futr = round(0.05/(1/fps));
    end
    if nargin < 6
        min_trans = round(0.1/(1/fps))-1;
    end
    fprintf('Obtaining features from dataset... \n');
    for m = 1:length(data) % For each csv file you uploaded.
        data{m} = data{m}*24/pix_cm; % With our zoom and resolution, our threshold was set according to the ~24 pixel/cm set-up.
        %% Obtain features
        clear fpd_norm cfp_pt_norm sn_pt_norm sn_pt_ang sn_disp pt_disp hpR_disp hpL_disp fpR_disp fpL_disp;
        cfp = [mean([data{m}(:,3),data{m}(:,5)],2),mean([data{m}(:,4),data{m}(:,6)],2)]; % Center of front paws x,y
        cfp_pt = [cfp(:,1) - data{m}(:,11), cfp(:,2) - data{m}(:,12)]; % Center of front paw to proximal tail x,y
        sn_pt = [data{m}(:,1) - data{m}(:,11), data{m}(:,2) - data{m}(:,12)]; % Snout to proximal tail x,y
        for i = 1:length(data{m}) % Euclidean distance of x,y, since position means nothing
            fpd_norm(i) = norm(data{m}(i,3:4)-data{m}(i,5:6)); % Front paw R to L euclidean distance    
            cfp_pt_norm(i) = norm(cfp_pt(i,:)); % Center of front paws to proximal tail euclidean distance
            sn_pt_norm(i) = norm(sn_pt(i,:)); % Snout to proximal tail euclidean distance, i.e. body length
        end
        fpd_norm_smth{m} = movmean(fpd_norm,[smth_hstry,smth_futr]); % Reduce label noise
        sn_cfp_norm_smth{m} = movmean(sn_pt_norm-cfp_pt_norm,[smth_hstry,smth_futr]); % Reduce label noise
        sn_pt_norm_smth{m} = movmean(sn_pt_norm,[smth_hstry,smth_futr]); % Reduce label noise
        for k = 1:length(data{m})-1 % Velocity and angle over time
            b_3d = [sn_pt(k+1,:),0]; a_3d = [sn_pt(k,:),0]; c = cross(b_3d,a_3d);
            sn_pt_ang(k) = sign(c(3))*180/pi*atan2(norm(c),dot(sn_pt(k,:),sn_pt(k+1,:))); % Body angle, arctan between body  
            sn_disp(k) = norm(data{m}(k+1,1:2)-data{m}(k,1:2)); % Snout displacement over time
            pt_disp(k) = norm(data{m}(k+1,11:12)-data{m}(k,11:12)); % Proximal tail displacement over time
            hpR_disp(k) = norm(data{m}(k+1,7:8)-data{m}(k,7:8)); % Right hind paw displacement over time
            hpL_disp(k) = norm(data{m}(k+1,9:10)-data{m}(k,9:10)); % Left hind paw displacement over time
            fpR_disp(k) = norm(data{m}(k+1,3:4)-data{m}(k,3:4)); % Right forepaw displacement over time
            fpL_disp(k) = norm(data{m}(k+1,5:6)-data{m}(k,5:6)); % Left forepaw displacement over time
        end 
        sn_pt_ang_smth{m} = movmean(sn_pt_ang,[smth_hstry,smth_futr]); % Reduce label noise
        sn_disp_smth{m} = movmean(sn_disp,[smth_hstry,smth_futr]); % Reduce label noise
        pt_disp_smth{m} = movmean(pt_disp,[smth_hstry,smth_futr]); % Reduce label noise
        hpR_disp_smth{m} = movmean(hpR_disp,[smth_hstry,smth_futr]); % Reduce label noise
        hpL_disp_smth{m} = movmean(hpL_disp,[smth_hstry,smth_futr]); % Reduce label noise
        fpR_disp_smth{m} = movmean(fpR_disp,[smth_hstry,smth_futr]); % Reduce label noise
        fpL_disp_smth{m} = movmean(fpL_disp,[smth_hstry,smth_futr]); % Reduce label noise
        %% Classify action based on parameters specified above
        %%% Rest/Pause is defined as minimal change in body angle, and minimal displacement of all labels.
        idx_rest{m} = find(abs(sn_pt_ang_smth{m}) < 1 & sn_disp_smth{m} > 0 & sn_disp_smth{m} < 2 & pt_disp_smth{m} < 2 & ...
            fpR_disp_smth{m} < 2 & fpL_disp_smth{m} < 2 & sn_cfp_norm_smth{m}(:,2:end) > 0);
        %%% Rear is defined as either snout is gone, front paws are gone (reared up), or front paws coming together (rearing up).
        idx_rear{m} = find(sn_disp_smth{m} == 0 & abs(sn_pt_ang_smth{m}) < 1 | ... 
            fpR_disp_smth{m} == 0 & abs(sn_pt_ang_smth{m}) < 1 | ... 
            fpL_disp_smth{m} == 0 & abs(sn_pt_ang_smth{m}) < 1 | ...
            sn_disp_smth{m} > 2 & abs(fpR_disp_smth{m}-fpL_disp_smth{m}) < 5 & sn_cfp_norm_smth{m}(:,2:end) > 0 & ... 
            fpd_norm_smth{m}(:,2:end) < 20 & pt_disp_smth{m} >= 2 & fpR_disp_smth{m} >= 2 & fpL_disp_smth{m} >= 2);
        %%% Groom is defined as either snout invisible, front paws visible and close, and no change in angle or tail position
        %%%% snout is visible, front paws in/visible and far, hands are in front of snout, and no change in angle or tail position
        %%%% or snout visibly moving, both front paws visible, one of which moves, and no change in angle or tail position
        idx_groom{m} = find(sn_disp_smth{m} == 0 & fpd_norm_smth{m}(:,2:end) < 10 & fpR_disp_smth{m} >= 2 & abs(sn_pt_ang_smth{m}) < 1 & ...
            pt_disp_smth{m} < 2 | sn_disp_smth{m} == 0 & fpd_norm_smth{m}(:,2:end) < 10 & fpL_disp_smth{m} >= 2 & ...
            abs(sn_pt_ang_smth{m}) < 1 & pt_disp_smth{m} < 2 | ...
            sn_disp_smth{m} > 0 & fpR_disp_smth{m} >= 2 & abs(sn_pt_ang_smth{m}) < 1 & hpR_disp_smth{m} < 2 & hpL_disp_smth{m} < 2 & pt_disp_smth{m} < 2 | ...
            sn_disp_smth{m} > 0 & fpL_disp_smth{m} >= 2 & abs(sn_pt_ang_smth{m}) < 1 & hpR_disp_smth{m} < 2 & hpL_disp_smth{m} < 2 & pt_disp_smth{m} < 2 | ...
            sn_disp_smth{m} > 0 & abs(sn_pt_ang_smth{m}) < 1 & sn_cfp_norm_smth{m}(:,2:end) < 0 & pt_disp_smth{m} < 2 | ...
            sn_disp_smth{m} >= 2 & fpR_disp_smth{m} > 0 & fpL_disp_smth{m} > 0 & abs(sn_pt_ang_smth{m}) < 1 & pt_disp_smth{m} < 2);
        %%% Orientation of body is defined as non-locomoting body angle change
        idx_ori_right{m} = find(pt_disp_smth{m} < 2 & sn_pt_ang_smth{m} >= 1);
        idx_ori_left{m} = find(pt_disp_smth{m} < 2 & sn_pt_ang_smth{m} <= -1);
        %%% Locomotion is defined as body displacement with either front paw
        idx_loc{m} = find(pt_disp_smth{m} >= 2 & fpR_disp_smth{m} >= 2 | pt_disp_smth{m} >= 2 & fpL_disp_smth{m} >= 2 | ...
            pt_disp_smth{m} >= 2 & hpR_disp_smth{m} >= 2 | pt_disp_smth{m} >= 2 & hpL_disp_smth{m} >= 2);
        %%% Head movement is defined as only snout movement with stationary body, and no change in angle or tail position
        idx_headmov{m} = find(abs(sn_pt_ang_smth{m}) < 1 & sn_disp_smth{m} >= 2 & hpR_disp_smth{m} < 2 & hpL_disp_smth{m} < 2);  
        %% Low pass filter of unreasonable fast state changes
        MC = []; MC(:) = 0; MC(idx_rest{m},1) = 1; MC(idx_loc{m},1) = 5; MC(idx_ori_left{m},1) = 6; MC(idx_ori_right{m},1) = 7; 
        MC(idx_rear{m},1) = 2; MC(idx_groom{m},1) = 3; MC(idx_headmov{m},1) = 4; 
        stch = find(diff(MC)~=0); % Remove short transitions
        for i = 1:length(stch)-1
            if stch(i+1) - stch(i) <= min_trans
               MC(stch(i):stch(i+1)) = MC(stch(i));
            end
        end
        %% Reassign sandwiched unidentified behavior to prior and posterior, if the same 
        for k = 1:length(stch)-2
            if MC(stch(k+1)) == 0 && MC(stch(k)) == MC(stch(k+2))
               MC(stch(k):stch(k+2)) = MC(stch(k));
            end
        end
        g_num{m} = MC; % Store behavior per .csv
        %% Pull out transitional probability
        perc_unk{m} = length(find(g_num{m}==0))/(length(data{m})-1); % Undefined percentage
        binEdges = 1:length(MC); [~,states] = histc(nonzeros(MC),binEdges); % Histc just asks which bin each pose belongs to
        binTicks = binEdges(1:end-1); % Don't need the last bin edge
        ActMat = zeros(length(unique(MC))-1); fromBinNos = states(1:end-1); nextBinNos = states(2:end);
        for i = 1:length(fromBinNos)
            ActMat(fromBinNos(i),nextBinNos(i)) = ActMat(fromBinNos(i),nextBinNos(i)) + 1;
        end
        for j = 1:length(unique(MC))-1
            trans_p{m}(j,:) = ActMat(j,:)/sum(ActMat(j,:));
        end
    end

return

